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The structure of liquid water at ambient conditions is studied in ab initio molecular dynamics 
simulations in the NVE ensemble using van der Waals (vdW) density- functional theory, i.e. using 
the new exchange-correlation functionals optPBE-vdW and vdW-DF2, where the latter has softer 
non-local correlation terms. Inclusion of the more isotropic vdW interactions counteracts highly 
directional hydrogen-bonds, which are enhanced by standard functionals. This brings about a 
softening of the microscopic structure of water, as seen from the broadening of angular distribution 
functions and, in particular, from the much lower and broader first peak in the oxygen-oxygen pair- 
correlation function (PCF) and loss of structure in the outer hydration shells. Inclusion of vdW 
interactions is shown to shift the balance of resulting structures from open tetrahedral to more close- 
packed. The resulting O-O PCF shows some resemblance with experiment for high-density water 
(A. K. Soper and M. A. Ricci, Phys. Rev. Lett., 84:2881, 2000), but not directly with experiment 
for ambient water. Considering the accuracy of the new functionals for interaction energies we 
investigate whether the simulation protocol could cause the deviation. An O-O PCF consisting of a 
linear combination of 70% from vdW-DF2 and 30% from low-density liquid water, as extrapolated 
from experiments, reproduces near-quantitatively the experimental 0-0 PCF for ambient water. 
This suggests the possibility that the new functionals may be reliable and that instead larger-scale 
simulations in the NPT ensemble, where the density is allowed to fluctuate in accordance with 
proposals for supercooled water, could resolve the apparent discrepancy with the measured PCF. 



I. INTRODUCTION 

Liquid water plays a crucial role in all biological and 
numerous chemical processes, which has provided the in- 
centive for many detailed experimental and theoretical 
studies probing both structural and dynamical proper- 
ties of the fluid. However, the microscopic structure 
of ambient liquid water is still a matter of current de- 
bate li"— In particular two classes of models are currently 
being considered, where the traditional model of water is 
based on a continuous distribution of distorted tetrahe- 
dral structures. This is typical of what most molecular 
dynamics simulations currently give. However, most of 
these simulations give over-structured 0-0 and 0-H pair- 
correlation functions (PCFs) and show discrepancies in 
comparison to x-ray and neutron scattering experimen- 
tal datat^ii^ It is, however, possible to generate a more 
distorted tetrahedral structure model that is consistent 
with the diffraction data, but equivalent agreement is 
seen also for alternative asymmetrical and mixture mod- 
els illustrating that diffraction data do not discriminate 
between differently hydrogen-bonded (H-bonded) struc- 
ture models 

Based on recent findings correlating x-ray emission 



spectroscopy (XES) with x-ray Raman scattering (XRS) 
and small angle x-ray scattering (SAXS) data^i^ a model 
has been suggested where a division into contributions 
from two classes of local instantaneous H-bonded struc- 
tures is driven by incommensurate requirements for min- 
imizing enthalpy and maximizing entropy; in particular 
the XES data show two well-separated peaks which in- 
terconvert but do not broaden with changes in temper- 
ature f^^'i^i^'S^ In the proposed picture the dominating 
class at ambient temperatures consists of a continuum 
of structures with some resemblance to high-pressure 
water fiS but with a further expanded first shell (more 
distorted H-bonds) and more disorder in the 2nd shell; 
this was based on the temperature dependent shift of 
the dominating peak in the XES spectra indicating more 
thermal distortion and disorder with increasing tempera- 
ture. The second class corresponds to fluctuations where 
regions of strongly tetrahedral structures (similar to low- 
density water) appear in different sizes and shapes as the 
molecules attempt to form enthalpically favored tetrahe- 
dral H-bond structures, resulting in mean size interpreted 
from the SAXS data as ^lum.^ but naturally many sizes 
and shapes would appear. It should be emphasized that 
since these are fluctuations no strict boundaries between 



the two classes should be expected. 

The attosecond (XRS, SAXS) to femtosecond (XES) 
time scales of the experimental probes are too fast for 
molecular motion to be followed and the experimental 
data thus correspond to a statistical sampling of instan- 
taneous, frozen local structures in the liquid; no exper- 
imental information on the time scale of such fluctua- 
tions is thus currently available.^- Besides being con- 
sistent with both neutron and x-ray diffraction/^ this 
picture was recently also shown to bring a consistency 
between x-ray diffraction and extended x-ray absorp- 
tion fine structure (EXAFS) data, requiring both disor- 
dered structures and a fraction of molecules with straight, 
strong hydrogen-bonds Other opinions, however, ex- 
ist regarding the interpretation of the new SAXS, XES 
and XRS data.-^'^^'i^^'Si-^ On the other hand, recent 
SAXS data extending into the supercooled regime and 
supported by theoretical simulations^!^ as well as re- 
cent high-quality x-ray diffraction data resolving shell 
structure out to 12 A even in ambient and hot water, 
but contributed by a minority species, "^^ provide support 
for the original interpretation in this debate. 

Most simple rare gas solids and liquids have a nearest- 
neighbor coordination of 12 whereas hexagonal ice, due 
to the directional H-bonds has a coordination of only 4. 
The latter leads to large open volumes in the ice lat- 
tice and a resulting low density. The dispersion, or van 
der Waals (vdW) force, in condensed rare gases leads to 
non-directional, isotropic interactions and closer packing. 
Similarly, the inclusion of vdW interactions in ab initio 
simulations of water may counteract the directional in- 
teractions and lead to better agreement with, e.g., ex- 
perimental PCFs. Here it should be understood that, 
while this could be regarded as a minimum requirement 
of a water model, it is by no means sufhcient for a com- 
plete description. Interestingly, it has been argued on 
thermodynamic grounds that over a large range of the 
liquid-vapor coexistence line the averaged water interac- 
tion potential should resemble that of liquid Argon^^, i.e. 
not be determined by directional H-bonding. 

Water shows many anomalies in its thermodynamic 
properties, such as compressibility, density variation and 
heat capacityi^"— In attempts to explain this, direc- 
tional H-bonds and more isotropic vdW forces are key 
concepts. While vdW forces are well defined as results of 
non-local electronic correlations, there is no unique way 
to characterize H-bonds in terms of topology or interac- 
tion strength. And yet "the H-bond governs the overall 
structure and the dynamics of water" 

One of the models to explain the enhanced anoma- 
lies in supercooled water is the liquid-liquid critical-point 
(LLCP) hypothesis,™"- with the most substantial role 
played by cooperative H-bond interactions among the wa- 
ter molecules4i The LLCP model explains the signifi- 
cant increase in density fluctuations upon supercooling 
water, which is evidenced by the anomalously increasing 
isothermal compressibility,4S. as resulting from attempts 
to locally form enthalpically favored open tetrahedrally 



coordinated H-bond regions. It furthermore connects the 
deeply supercooled liquid state of water to the polyamor- 
phism seen in ices, i.e. the distinct low-density and 
high-density amorphous ice phases (LDA/HDA). A high- 
density liquid (HDL) phase transforms to an ordered low- 
density phase (LDL) in the deeply supercooled region 
through a first-order phase transition at high pressures 
above the LLCP and through a continuous smooth tran- 
sition upon crossing the Widom line at pressures below 
the criticaLiS"^"-^ There are differences in their respec- 
tive local structures; in pure HDL the local tetrahedrally 
coordinated H-bond structure is perturbed by a partially 
collapsed second coordination shell, while in the LDL a 
more open and locally " bulk-ice-like" H-bond network is 
realized.iii^'^ 

The combined XES, XAS and SAXS resuhs described 
above,— which indicate nanoscale density and structural 
fluctuations, can be easily interpreted as reflections of 
this "competition" between the two local forms, HDL 
(maximizing entropy) and LDL (minimizing enthalpy) 
and thus viewed as extending an established picture of su- 
percooled water into the ambient regime. Whether HDL 
and LDL can exist as pure phases, accompanied by a 
liquid- liquid phase transition and a critical point, is still 
unresolved and alternative models, e.g., singularity-free 
(SF)^ii^ critical-point-free (CPF)43 and stability hmit 
(SL) conjecture^ scenarios have been proposed, however 
still building on structural HDL/LDL fluctuations. 

In the quantitative characterization of water, computer 
simulations play a vital role. Empirical force fields are 
frequently applied but with questionable transferability, 
since force fields are parameterized against experimental 
data or against a by necessity limited set of quantum 
chemically computed structures. Furthermore, many- 
body interactions beyond pair-interactions are frequently 
not taken into account. 

These deficiencies are eliminated in Car-Parrinello^ 
(CP) and Born-Oppenheimer (BO) molecular dynamics 
(MD), collectively known as ab initio (AI) MD. In AIMD, 
the forces are calculated using a first-principles electronic 
structure method, typically based on density functional 
theory (DFT). BOMD, used in the present study, min- 
imizes the Kohn-Sham energy functional at each time 
step, keeping the nuclear positions frozen. In nearly all 
force field and AIMD simulations of water at ambient 
conditions there seems to be a strong driving force to 
form highly directional H-bonds, leading to tetrahedral 
structures that are in general over-structured in terms of 
the derived PCFs. One exception is the coarse-grained 
mW water model^^, which has two terms in the inter- 
action potential corresponding to anisotropic tetrahedral 
interactions and isotropic vdW interactions, respectively, 
and which gives a maximum peak height of 2.3 in the 0-0 
PCF at room temperature, in close agreement with recent 
analyses of experimental diffraction results.— 'i^'^'^"— 
This model was shown to feature fluctuations between 
tetrahedral and disordered species resulting in a liquid- 
liquid transition in the supercooled region. Empirical 
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force-field models which have over-structured PCFs in 
agreement with older determinations^!^ have however 
also been shown to exhibit liquid-liquid phase transitions 
in the supercooled regime, e.g., Refs. I40lf59ll60l and IsTI. in- 
dicating that the PCFs are not decisive for general trends 
in the thermodynamic behavior in water simulations. 

Until recently, AIMD simulations of water have almost 
exclusively been performed with the BLYP^^ and PBE^'^ 
exchange-correlation (XC) functionals. However, these 
functionals are shown to significantly over-structure liq- 
uid water as seen from the maximum value and sharp- 
ness of the first peak in the oxygen-oxygen PCF com- 
pared to recent data and analyses.— '^^'2^'^"— AIMD 
simulations of water have furthermore been shown to de- 
pend on which functional is applied and to give differ- 
ent predictions for different XC functionals.^^ MD sim- 
ulations performed using functionals based on the gen- 
eralized gradient approximation (GGA) tend to over- 
structure liquid water and lead to diffusion constants 
two to three times too small compared to experiment; 
using hybrid functionals only marginally improves the 
results;2S In addition, it has been shown that PBE-based 
AIMD simulations lead to a melting point of ice at 417 
K and therefore simulations at ambient conditions with 
this functional will describe a deeply supercooled state 
which is strongly over-structured with respect to real liq- 
uid water at ambient conditions. As we will show in the 
following, inclusion of the more isotropic vdW interaction 
balances the directional forces allowing a partial break- 
down of the H-bond network and a much less structured 
liquid. 



II. METHODS 
A. Role of van der Waals (vdW) forces 

Small water clusters have been studied using the PEE 
and BLYP XC functionals, which do not explicitly in- 
clude vdW interactions, and the results compared to high 
accuracy methods such as coupled cluster (CCSD(T)) 
and M0ller-Plesset (MP2). With PBE^^ near chemical 
accuracy for the strength of the H-bond for the water 
dimer is obtained while BLYP consistently underbinds 
small water clusters.-^ However, discrepancies arise and 
increase with the size of the water cluster for both PBE 
and BLYP. This has been ascribed to the lack of a de- 
scription of vdW forces 1^ One could thus argue that ob- 
taining the correct result for the water dimer is essential 
but no guarantee for a correct description since not all 
physical interactions relevant for larger clusters are sam- 
pled by the dimer. 

While it is well established that at low temperatures 
H-bonds give the major contributing factor to the dy- 
namics and structure of water, vdW interactions have 
also been suggested to be importan t ''"i^^ . In line with 
this, thermodynamic considerations have led to the sug- 
gestion that at higher temperatures the averaged water 



interaction potential should resemble that of liquid Ar- 
goni.22, The angular dependence of the H-bond is antici- 
pated to have a big impact on the PCF and self-diffusion 
coefficient.''^- If, for example, it is too difficult to bend a 
DFT H-bond, the diffusion coefficient should come out 
too small, which it does. Many other suggestions to ex- 
plain the too small diffusion coefficient exist however— 
but balancing the directional H-bond interactions with 
more isotropic vdW forces would intuitively contribute 
to softening the H-bond network and allow more efficient 
diffusion. Traditional local and semi-local DFT do not, 
however, contain non-local vdW interactions, e.g., BLYP 
being especially incapable of describing dispersion*?'^ In- 
fluences of vdW interactions have been investigated using 
MD based on empirical potentials, ^^'^''^ e.g., performed 
with a dispersion-corrected BLYP XC functional/^ or 
using empirically damped CgR^® corrections^^"— to de- 
scribe the vdW interactions. 

A way to introduce vdW forces in DFT from first prin- 
ciples is provided by the van der Waals density functional 
vdW-DF^ recently used for the first time in AIMD on 
liquid water The inclusion of vdW forces using the 
vdW-DF was shown to greatly improve water's equilib- 
rium density and diffusivity. However the vdW-DF MD 
also produces a collapsed second coordination shell giving 
rise to new structural problems that have been suggested 
to depend partially on the choice of exchange used in the 
vdW functional.*^ 

The vdW-DF method proposed by Dion et aZ,~ ac- 
counts for exchange by a functional that gives Hartree- 
Fock-like repulsion at relevant separations and that in- 
cludes non-local correlation, and thus vdW forces, by 
calculating the dielectric response in a plasmon-pole ap- 
proximation. It gives the correct stability trend for low- 
lying water hexamers*^ but returns a significant under- 
binding for most H-bonds.—^— The underbinding can 
be remedied by using an exchange functional that gives 
more binding*^ at typical H-bond separations f^^iS^i^ like 
the PW86,Si optPBE^ and C09*^^ exchange functionals. 
Recently Klimes et al^ proposed a new vdW density 
functional, optPBE-vdW, based on the original vdW- 
DF functional.*" This scheme shows promise in the de- 
scription of dispersion and H-bonded systems, as it re- 
duces the underbinding given by the vdW-DF down to 
chemical accuracy while preserving the correct hexamer 
trends. However, this improved behavior is obtained at 
the cost of poorer performance on the binding energy 
of small moleculesi^ Very recently a second version of 
the vdW-DF, called vdW-DF2, was suggested,^- using a 
new non-local correlation functional along with a slightly 
refitted version of the PW86, called PW86R*^ as an ap- 
propriate exchange functional. Both optPBE-vdW and 
vdW-DF2 give chemical accuracy for the water dimer, 
albeit with slightly different balance between non-local- 
correlation and exchange contributions. In the present 
study we therefore wish to investigate the microscopic 
structure of liquid water by performing AIMD using both 
the new optPBE-vdW and vdW-DF2 XC functionals to 
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also investigate the importance of the balance between 
correlation and exchange in liquid water AIMD simula- 
tions. 



B. The optPBE-vdW and vdW-DF2 
exchange-correlation functionals 

In general a vdW-DF functional takes the form 

where E^'^^ is an exchange functional using the general- 
ized gradient approximation (GGA), E"^^^ accounts for 
the local correlation energy by using the local density 
approximation (LDA). LDA is chosen to avoid double 
counting of correlation. The non-local correlation en- 
ergy describing the vdW interaction is given by the six- 
dimensional integral^ 

Ef = \j j n{v)cjy{v,v')n{v')dvdv', (2) 

where 0(r, r') is the interaction kernel and depends on the 
density and its gradient. The non-local term is calculated 
as suggested in Ref. lo^. In the original vdW-DF from 
Dion et al. the exchange functional from revPBE-^'^ is 
utilized. 

The optPBE-vdW functional is constructed like vdW- 
DF— but uses an alternative exchange functional. The 
latter takes the same form as both the PBE and RPBE 
exchange, but the parameters of the exchange enhance- 
ment factor are optimized against the S22 database.*^ 
The S22 database^ is a set of 22 weakly interacting 
dimers, mostly of biological importance, including the 
water dimer. 

The vdW-DF23i has the form of Eq. (l)and uses the 
PW86 exchange f25 which is argued in Ref. [Sa to give the 
most consistent agreement with Hartree-Fock (HF) ex- 
act exchange, and with no spurious exchange binding. 
Furthermore, a new approximation for E"' is used to cal- 
culate the value of the interaction kernel in Eq. (2))^ 
This new functional has been shown to give very accu- 
rate results for the water dimer as compared to bench- 
mark CCSD(T) calculations2ii2S and to compare closely 
to the S22 benchmarki22 



eV per electron is used in order to determine the forces 
adequately. 

All internal bond lengths are kept fixed at 0.9572 A 
(an MP2 optimized gas phase geometry obtained from 
the G2-database)^™ but angles are allowed to vary (i.e. 
bending vibrations are included); eliminating the high- 
frequency OH-stretch allows longer time steps in the sim- 
ulations albeit introducing some uncertainty^Si which, 
however, is not relevant for the large differences observed 
in our simulations between the PBE on the one hand 
and the vdW functionals on the other since all simula- 
tions have this constraint imposed. In the initial config- 
uration 64 water molecules are placed in a simple cubic 
lattice with random orientations in a cubic periodic box 
with side lengths 12.42 A, to reproduce a water density 
of 1 g/cm^. The geometry is then optimized to obtain 
a configuration at zero Kelvin (using PBE), from which 
the MD is started giving the atoms random velocities 
according to a Maxwell-Boltzmann velocity distribution 
corresponding to two times 300K, keeping the center of 
mass of the box stationary. Approximately half of the 
kinetic energy converts to potential energy thus giving 
an average temperature around 300K. An initial equili- 
bration of 10 ps using the PBE XC functional is per- 
formed followed by 2.5 ps vdW equilibration of the simu- 
lations using optPBE-vdW and vdW-DF2. For all meth- 
ods equilibration was followed by production runs for 10 
ps which is the minimum time reported necessary due to 
the slow diffusion of water Using 64 water molecules 
has been shown to be adequate to remove the most signif- 
icant problems concerning finite size effects^S^ and is fea- 
sible within the current computational capabilities. The 
Verlet algorithm is employed using a time step of 2 fs 
in the NVE ensemble. Using this type of ensemble the 
temperature is allowed to fluctuate and the average tem- 
perature of the PBE, vdW-DF2 and optPBE-vdW sim- 
ulations were 299K, 283K and 276K, respectively. The 
same computational setup has been used for the PBE and 
vdW density functional MD simulations in order to allow 
direct comparison of the different models. Since simula- 
tions with PBE at ambient conditions describe a deeply 
supercooled state relative to its melting point at 417 K— 
the PBE simulations are only performed here to provide 
a reference for the effects of including vdW interactions 
through the optPBE-vdW and vdW-DF2 functionals. 



C. Computational protocol 

Ab initio molecular dynamics simulations are per- 
formed in the NVE ensemble with optPBE-vdW, vdW- 
DF2, and PBE, using the grid-based real-space projector 
augmented wave GPAW codcj ^^i^^ A wave function grid 
spacing of 0.18 A and Fermi smearing with a width of 
0.01 eV have been used. The grid spacing has been de- 
termined by comparing DFT calculations of water hex- 
amers with CCSD(T) results. In the electronic structure 
calculations a strict energy convergence criterion of 10~^ 



III. RESULTS 

A. Water dimer 

Before discussing the MD results we compare the func- 
tionals for a simpler but still relevant system: the water 
dimer. Fig. [T]a) illustrates the potential energy curve for 
the water dimer calculated using PBE, vdW-DF, vdW- 
DF2 and optPBE-vdW in comparison with the bench- 
mark CCSD(T) curve from Ref. lOa. Fig. [T]a) shows that 
the vdW functionals are capable of describing this basic 
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FIG. 1: a) The water dimer potential energy curves calcu- 
lated with DFT using the XC functionals PBE, vdW-DF, 
vdW-DF2 and optPBE-vdW, respectively, are compared to 
CCSD(T)/CBS wave function results;^ b) The distance de- 
pendence of the non-local contribution (Eq. (2)) to the inter- 
action energy of the water dimer for the XC functionals vdW- 
DF2 and optPBE-vdW, shows that, when they give similar 
potential energies (Fig. la)) they do so for different reasons; 
the optPBE-vdW gets more binding from a stronger vdW at- 
traction but vdW-DF2 gets more net attraction from a less 
repulsive exchange. 



constituent of liquid water extremely accurately, however 
for different reasons. The non-local contribution (E"') to 
the dimer binding from the two functionals is plotted in 
Fig. [T]b). The non-local part of the optPBE-vdW func- 
tional, which is based on the older approximation, is more 
attractive as mentioned in Ref. 91. Since less attraction 
stems from the non-local interaction in the vdW-DF2, 
while the total energy for the dimer is almost identical 
to that of optPBE-vdW, the remaining part of the in- 
teraction energy must give a larger contribution for the 
vdW-DF2 than for optPBE-vdW. The remaining part 
of the interaction energy includes electrostatic interac- 
tion, electronic correlation, and more or less repulsive 
exchange. Since electrostatic interactions only depend 
on separation, and local correlation is treated identically 
with the LDA correlation in both cases, this difference 
has to come from the different choices for the exchange. 
The PW86 exchange in vdW-DF2 is hence less repulsive 
than the optPBE exchange in optPBE-vdW; a possible 
cause of the reported collapsed second-shell structure was 
in Ref. [Slj suggested to be that the non-local parameter- 
ization of exchange used in vdW-DF and optPBE-vdW 
may be too attractive when used in MD. This is, however, 
not the case, as seen from the pair-correlation functions 
(PCFs), to be discussed next. 



B. Pair-correlation functions 

Fig. [2] a) illustrates that AIMD simulations of liquid 
water using vdW-DF2 and optPBE-vdW give very sim- 
ilar 0-0 PCFs which are, however, very different from 
the 0-0 PCF from PBE and furthermore from those de- 
rived from experiment using either Empirical Potential 



Structure Refinement (EPSR)^^ or Reverse Monte Carlo 
(RMC)i04 ^i^g structure factorJJ^ In the simula- 

tions both functionals result in the same characteristics 
as reported in Ref. Hi], including a lower first peak shifted 
to larger 0-0 separation than for normal GGAs as well 
as for experiment on ambient water. The second coor- 
dination shell at 4.5 A is also completely smeared out 
where correlations from the region 4-5 A have instead 
moved into the region 3.3-3.7 A. The non-local correla- 
tion differences in the functionals do not, however, re- 
sult in significantly different 0-0 PCFs, but we note a 
somewhat higher (2.5) first peak for vdW-DF2 compared 
to optPBE-vdW (2.3). Since the latter gives a slightly 
stronger non-local contribution we take this as indication 
that it is indeed the vdW contribution that so strongly 
affects the first shell structure in the simulations. 

In contrast, the very recent vdW-DF MD simulation 
showed that by changing the exchange in vdW-DF from 
revPBE to PBE, the second shell structure again be- 
came well defined.— However, the exchange functionals 
of revPBE and PBE are quite different, making an expla- 
nation in terms of the exchange less likely; the potential 
energy curve of the dimer is furthermore not reproduced 
very well using the PBE exchange with LDA and non- 
local correlation suggesting that substituting revPBE by 
PBE for the exchange does not lead to consistent im- 
provement in the description. 
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FIG. 2: a) Oxygen-oxygen PCFs (goo) obtained from exper- 
imental data using EPSR^'' and RMC^^ in comparison with 
PCFs obtained by DFT MD simulations using PBE, optPBE- 
vdW and vdW-DF2. b) Oxygen-hydrogen PCFs (gou) ob- 
tained from experimental data using EPSR^ and RMOi^ m 
comparison with PCFs from PBE, optPBE-vdW and vdW- 
DF2. 



Compared to the experimentally derived 0-0 PCFs it 
is clear that the PCF obtained from PBE is severely over- 
structured while the simulations including vdW forces 
have resulted in a significantly less structured PCF than 
what is experimentally observed for ambient liquid wa- 
ter. Clearly neither simulation model gives direct agree- 
ment with the experimental 0-0 PCF even though, in 
the case of the vdW functionals, small water clusters are 
described very accurately. We will address this aspect in 
the discussion section below. 

Some discrepancy in the 0-H PCF for the vdW XC 
functionals compared to experiment is seen in Fig. [2]b), 
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which we shall now consider. The 0-0 correlations can 
be obtained from a Fourier transform of the x-ray diffrac- 
tion data, if a large enough k-range has been measured 
and the data can be properly normalized; x-ray scatter- 
ing is strongly dominated by the electron-rich oxygens. 
Neutron diffraction data, on the other hand, contain si- 
multaneous information on the H-H, 0-H and to some ex- 
tent, the 0-0 PCFs, making a direct Fourier transform 
to extract a specific PCF inapplicable. Various fitting 
schemes of structure models to the experimental struc- 
ture factors have therefore been developed, and we show 
two such fits to the same experimental data using the 
and RMC methodji^ respectively. 
There is a significant difference in the first peak posi- 
tion in the 0-H PCF between the EPSR and RMC fits 
compared to what is found for the 0-0 PCF. This can 
be understood from the relatively lower sensitivity of the 
neutron data to specifically the 0-H correlation in com- 
parison to the sensitivity of x-ray data to the 0-0 corre- 
lationji^ The EPSR technique uses the assumed reference 
pair-potential to provide structural aspects not included 
in the experimental data /^i^^ while structural aspects 
not determined by the experimental data, or imposed 
constraints, will in the RMC technique simply result in 
a phase-space weighted sampling of structures consistent 
with the experimental structure factors; combining the 
two methods thus gives additional information on the un- 
certainties and assumptions in the resulting fits. It is in- 
teresting to note that the RMC method gives a shift in 
the first peak of the 0-H correlation out to nearly 2 A,— 
which agrees well with the vdW MD simulations pre- 
sented here, while the EPSR solution is closer in position 
to the PBE, likely reflecting the SPC/E starting force- 
field in the EPSR fitting procedure. Note, that both the 
RMC and EPSR fits reproduce the experimental scatter- 
ing data equally well, implying that the position of the 
first intermolecular OH correlation is not strictly deter- 
mined by the data, which leaves an uncertainty in the 
diffraction-derived 0-H PCFJ^iiSe j.^^ gj.g^ pg^^j^ 
PBE 0-H PCF is clearly too high and the first minimum 
at 2.5 A too low, however, while all three simulations 
exaggerate the height of the second peak at 3.2 - 3.6 A; 
this can, however, be expec ted t o be reduced by including 
quantum effects, e.g., Ref. llOTl 



C. Angular distribution functions and 
hydrogen-bonding analysis 

The van der Waals functionals provide a smoother an- 
gular structure with less tetrahedral bonding as demon- 
strated by the angular distribution functions and the av- 
erage number of H-bonds per water molecule; here we 
use the cone criterion from Ref. d as a geometric H- 
bond definition: roo < r^o " 0.00044(5|oo- This de- 
fines a cone around each H-bond-donating OH group, 
where Tqq^ = 3.3 A is the maximum 00 distance at 
zero angle 5hoO: where (5hoo is the H-0- • -O angle quan- 



tifying the angular distortion of the H-bond. Table |T] 
shows the H-bond statistics for PBE, optPBE-vdW and 
vdW-DF2. PBE is seen to prefer a tetrahedral H-bond 
coordination with a majority of the molecules having 4 
H-bonds. Including non-local correlation has a large ef- 
fect where, for both optPBE-vdW and vdW-DF2, the 
H-bond distribution shifts from a majority with four H- 
bonds to instead a predominance of species with two or 
three. Comparing the two vdW functionals we observe 
that the optPBE-vdW has a slightly larger amount of wa- 
ter molecules having two or three H-bonds compared to 
vdW-DF2 which we ascribe to the relatively more repul- 
sive exchange and stronger non-local contribution in the 
former. The vdW-DF2 with its relatively weaker vdW 
interaction shows slightly higher preference to forming H- 
bonds. This analysis suggests that there is a competition 
between isotropic vdW forces and directional H-bonds, 
resulting in fewer or more H-bonds per water molecule 
depending on the applied approximations; however, be- 
tween the vdW models the average number of H-bonds 
varies only weakly despite differences in vdW strength. 
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FIG. 3: Molecular angular distributions in liquid water ac- 
cording to MD simulations using DFT with the indicated 
functional, a) The angular distribution functions of the O- 
H- . .0 angle, b) the H-O- • -O (first peak) and 6 =<H- ■ -O-H 
(second peak) angles obtained using the XC functionals PBE 
(yellow), vdW-DF2 (red) and optPBE-vdW (green). When 
including vdW interactions a softening of the structure is seen 
from the broader distribution of angles. 

We note in particular the low number of double-donor, 
double- acceptor tetrahedral molecules according to the 
cone criterion-^ for the two vdW models. In fact, the large 
number of broken H-bonds in the vdW simulations sug- 
gests that these models are in closer agreement with pre- 
dictions from x-ray spectroscopie a^'^'^'^i^^'^^" — compared 
to most other AIMD models and future calculated x-ray 
spectra based on optPBE-vdW and vdW-DF2 structures 
may provide an interesting opportunity to obtain further 
insight regarding the interpretation of these spectra. 
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No. H-bonds \ Method 


PBE 


optPBE-vdW 


vdW-DF2 


1 


2 


10 


8 


2 


12 


29 


27 


3 


31 


37 


38 


4 


52 


22 


25 


5 


3 


1 


1 



as 



TABLE I: Percentage distribution of hydrogen bonds per wa- 
ter molecule calculated using the cone criterion from Ref. Q. 
PBE results are shown to favor more H-bonds compared to ei- 
ther vdW functional, which both allow for a larger number of 
molecules to break the tetrahedral structure with four bonds. 



The angular distribution functions (ADFs) of the El- 
bonds are shown in Fig. [3] The ADFs of H-bond ac- 
ceptor and donor give information on the orientational 
flexibility of the water molecules. In the ADFs only the 
angles between a central molecule and the molecules of 
the first solvation shell are considered by using a cutoff 
distance corresponding to the first minimum in the PBE 
0-0 PCF; this distance was applied also to the vdW 
MDs where the second shell is smeared out and no min- 
imum is visible. Fig. |3] a) displays the distributions of 
donor angles a =<0-H- • -O for the various simulations. 
The first peak in Fig. [3]b) is (3, the deviation of the O- 
H- • -O bond from being linear, which gives information 
on the flexibility of donor H-bonds, while the second peak 
is the acceptor angle 6 =<H- • -O-H. The distribution of 
angles has been found to depend on the choice of wa- 
ter modelii^ The picture of a competition between non- 
directional vdW interactions and directed H-bonds seems 
to be supported by the ADFs as illustrated by the fact 
that the model without vdW forces (PBE) has no incen- 
tive to deviate from a structure of strong H-bonds, thus 
resulting in a relatively straight H-bond angle. When 
including vdW forces the H-bonds become significantly 
more bent. In general a softening of the structure is seen 
from the broader ADFs obtained in case of the vdW-DFs. 



D. Tetrahedrality and asphericity 

Two useful measures of the local coordination of 
molecules in water are the tetrahedralityi22iii2S and as- 
phericityii^ parameters. The former quantifies the de- 
gree of tetrahedrality in the nearest neighbor 0-0-0 an- 
gles and is defined as 



1=1 j>z 



COS Oioj + i 



(3) 



where 9ioj is the angle formed by two neighboring oxygen 
atoms i and j and the central molecule 0. Only the four 
nearest neighbors are taken into account which makes Q a 
very local measure. Perfect hexagonal ice gives (5 = 1 for 
all molecules while the ensemble average over an ideal gas 
gives < Q >= 0ii2£ The asphericity parameter is defined 



77 : 



367rT/2 



(4) 



where A and V are the area and volume of the Voronoi 
polyhedron of the molecule in question. Contrary to Q, 
rj is sensitive also to interstitial molecules outside the 
first shell and to the second coordination shell since these 
add surfaces to the Voronoi polyhedron, making it more 
spherical. The two relevant limits for water are that of 
hexagonal ice, which gives 77 — 2.25, and that of a perfect 
sphere which gives 77 = 1; larger disorder in the local 
coordination thus gives smaller values of 77. 
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FIG. 4: a) Distributions of the tetrahedrality parameter Q. 
vdW interactions lead to significantly lower average tetrahe- 
drality and a strong low-Q peak from interstitial molecules 
around Q = 0.5. b) Distributions of the asphericity parame- 
ter rj. A large effect of vdW interactions is seen with a shift 
towards more spherical (less ice-like) values. 

As Fig. |4] shows, the inclusion of the vdW interac- 
tion not surprisingly has a dramatic effect on both the 
tetrahedrality and asphericity distributions. The PBE 
simulation displays a strong peak at Q = 0.8, signifying 
a dominance of locally tetrahedral 0-0-0 angles, while 
both vdW simulations show an attenuation and shift of 
the high-(5 peak to lower tetrahedrality along with the 
appearance of a strong low-Q peak associated with in- 
terstitial molecules at non-tetrahedral positions between 
the first and second coordination shells. Out of the two 
vdW models, optPBE-vdW is seen to be somewhat less 
tetrahedral, consistent with their differences in H-bond 
statistics and PCFs discussed above. This is clearly illus- 
trated by the average tetrahedrality which is 0.692, 0.602 
and 0.583 for PBE, vdW-DF2 and optPBE-vdW, respec- 
tively. In comparison the average tetrahedrality has been 
estimated to be 0.576 using the EPSR methodjiii note, 
however, that the tetrahedrality parameter is experimen- 
tally rather uncertain, i. e. the same diffraction data have 
been shown to support tetrahedrality values ranging from 
0.488 to 0.603.1^ 

An even larger difference is seen in the asphericity dis- 
tributions; the two vdW models show sharper peaks cen- 
tered at lower asphericity values compared to PBE. This 
directly reveals the large disorder in second-shell corre- 
lations in the vdW models, resulting from the tendency 
to form more isotropic local structures when vdW forces 
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are included. Similarly to the comparison between the 
PCFs of the two vdW models discussed above, it can 
be seen here that despite non-local differences between 
vdW-DF2 and optPBE-vdW their respective liquid wa- 
ter structures turn out to be rather comparable in terms 
of both first- and second-shell correlations. The average 
asphcricity is f .681, 1.552 and 1.552 for PBE, vdW-DF2 
and optPBE-vdW, respectively. 



IV. DISCUSSION 

Considering the accuracy of the present versions of 
non-local correlation functionals, as calibrated against 
benchmark CCSD(T) and MP2 calculations for water 
dimer, water hexamers,— the S22 database,— we will 
here explore the possibility that the interactions between 
molecules in the simulation box are given sufficiently ac- 
curately by the functionals and that the resulting dis- 
crepancy between simulated and observed 0-0 PCF is 
rather due to limitations and constraints in the simula- 
tion protocol. 

Comparison of the results from the simulations using 
PBE with those including the vdW interactions shows a 
strong shift in the balance between directional H-bonding 
and more isotropic interactions; the former leads to tetra- 
hedral H-bond coordination and low density while the 
latter favors a more close-packed ordering and higher 
density, as evidenced by the loss of distinction between 
first and second coordination shells and the reduced num- 
ber of H-bonds. The simulations have in all cases been 
performed with internal OH distances fixed to the gas 
phase value; eliminating the high-frequency OH stretch 
allows longer time-steps to be used in the AIMD, but 
not allowing the internal OH distance to vary according 
to H-bond situation has been shown to lead to somewhat 
less structured PCFs in earlier worki^Si However, since 
the simulations with PBE, optPBE-vdW and vdW-DF2 
were all run with the same constraint in terms of internal 
OH distance this cannot explain the large effects on the 
0-0 PCF from including the vdW non-local correlation. 

We note that recent, high-precision x-ray diffraction 
measurements^ of ambient (25°C) and hot (66°C) wa- 
ter resolve shell-structure out to ~12 A in agreement 
with conclusions from SAXSi^; shell structure out to the 
fifth neighbor distance has been resolved before but only 
for supercooled water. ^^'^^^ Based on analysis of large- 
scale simulations with the TIP4P /2005 force- fieldii^ the 
shell structure could be assigned as due to an instanta- 
neous LDL-like minority species;22. The observed spatial 
extent of the correlation (12 A) is similar to the size of 
the present simulation box (12.42 A) making it unlikely 
that the simulation box is sufficiently extended to sup- 
port such experimentally observed instantaneous struc- 
tures. 

We compare the PCFs from the simulations performed 
using the vdW-DF2 and optPBE-vdW XC functionals 
(Fig. [5] a)), respectively, with the results of a neutron 



diffraction study*^ where LDL and HDL 0-0 PCFs were 
extrapolated from data at different pressures; the result- 
ing PCFs are shown in Fig. [5]b). The EPSR derived 
HDL PCF is rather similar to the PCF obtained using a 
Fourier transform of x-ray diffraction data at high pres- 
sures^^^ and furthermore seen to be very similar in terms 
of the second- and third-shell structure to that derived 
from vdW-DF2 and optPBE-vdW MD simulations; the 
effect of increasing pressure on the 0-0 PCF is that the 
4.5 A correlation disappears and moves to the 3.3 - 3.7 
A region and the third shell is shifted down to 6 AJ^ 
The 0-0 PCFs obtained using the vdW functionals sim- 
ilarly show a lack of well defined structure at 4.5 A, an 
increase in correlations at 3.3 - 3.7 A and show a shift 
towards shorter separations in comparison to PBE of the 
correlation at 6 - 6.5 A, as is seen from Fig. [5] Both are 
clear indications towards HDL water. However, in con- 
trast to the high pressure PCFs, a well defined peak at 
3.5 A is not present in the vdW MD simulations, but only 
an increase in correlations, and the first peak position is 
shifted outwards, which is not observed for pressurized 
water. 

Assuming that the AIMD simulations with non-local 
correlation and more isotropic interactions have led to 
a more compact, HDL-likc structure, it could be argued 
that a well-defined peak at 3.5 A should not be expected 
since, as deduced from XES spectra at different temper- 
aturesi^ii^, HDL-like water at ambient conditions should 
be thermally excited with a more expanded first shell 
and therefore further disordered in comparison to HDL 
water obtained under pressure. In particular, entropy 
effects due to thermal excitations leading to higher dis- 
order can be expected to create a structure where both 
the first shell and, in particular, the collapsed second 
shell are distributed over a range of distances, leading to 
molecules in what is often denoted interstitial positions 
and with the first 0-0 peak appearing at longer distance 
when not under pressure. In this respect a comparison 
with the amorphous high-density (HDA) and very high- 
density (VHDA) ices is of interest, where, for VHDA, 
the second shell moves inwards and a peak at 3.4 A 
develops while for HDA a peak is found at 3.7 A and 
the second peak broadens significantly. This indicates 
that various interstitial sites may be occupied making 
the high-density forms less well-definedJ^^-i^ It should 
be mentioned that a peak at ~3.7 A is present in the MD 
simulation performed by Wang et alS^ using the earlier 
vdW-DF ^'^ formulation of the functional. 

If we consider the proposed model of fluctuations be- 
tween HDL and LDL,i'2.ii2i it could well be that the 
vdW models under the present conditions only gener- 
ated HDL-like structures while without including vdW 
the resulting structure is clearly more LDL-like. Having 
two balancing interactions that favor opposite structural 
properties is a prerequisite for fluctuations; it is clear that 
by tuning either the importance of H-bonding or the vdW 
interaction the preference for either structure will be af- 
fected in the simulations. However, if the two proposed 
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structures of liquid water truly do coexist as endpoints of 
fluctuations in nanosized patches of different local den- 
sity, as suggested in Ref. [la, then an AIMD with only 64 
water molecules in a fixed volume may not be suitable 
to observe this behavior; a much larger box size and an 
NPT ensemble simulation allowing the box size to vary 
would be required. The relatively small simulation (12.42 
A box length) and short run time (10 ps) may only ob- 
serve a local structure of water which, in this picture, 
is either approximating LDL- or HDL-like. It should be 
noted that the simulations are run in the NVE ensemble 
with density fixed to correspond to ambient conditions 
which, under the assumption that ambient water is dom- 
inated by HDL, should furthermore favor an HDL-like 
structure over fluctuations towards LDL, if energetically 
allowed, as seems to be the case with vdW interactions 
included. 
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FIG. 5: a) Oxygen-oxygen PCFs (goo) obtained by MD sim- 
ulations from DFT with PBE, optPBE-vdW and vdW-DF2 
functionals. b) Experimental PCFs for high- and low-density 
water.— 



Exploring the hypothesis that the experimentally mea- 
sured 0-0 PCF in reality is the result of a spatial or tem- 
poral ave rag e over fluctuating structures as suggested in, 
e.g., Rei.M, and that the vdW-DF2 and optPBE-vdW 
functionals actually provide a sufficiently accurate inter- 
action potential, we will consider what additional con- 
tribution would be required to achieve agreement with 
the measured 0-0 PCF. The PCFs are, however, not 
directly measurable but derived from experimental data 
and we first need to discuss specifically the choice of O- 
O PCF for the comparison since different reference PCFs 
are used in the literature. 

X-ray and neutron diffraction data treated in conjunc- 
tion, either by the technique of empirical-potential struc- 
ture refinement (EPSR)^ or by reverse Monte Carlo 
(RMC) simulations r^ii^ as well as from directly Fourier 
transforming the latest high-quality x-ray diffraction 



data sets22i^i^ and the early data set of Narten and 
coworker&i22ii2^ all give a broad and slightly asymmet- 
ric first 0-0 peak with height 2.1-2.3 which is signifi- 
cantly lower than from standard MD simulations (height 
^3) and from previous analyses of either only neutron 
diffraction data using EPSR^ or analysis of the total x- 
ray scattering I(k) in terms of comparison to computed 
I(k) from MD simulations i^^i^ 

There were, however, problems with both the latter 
approache d^ '''^^1^^^ since neutron diffraction mainly mea- 
sures H-H and 0-H correlations and thus contain insuffi- 
cient information to modify the initial SPC/E force-field 
guess in EPSR to a solution that also describes the 0-0 
PCF, which is mainly determined by x-ray diffraction. 
The assumption by Hura et ^^ ,58,124 ^^^^ ^-^^^^ some ex- 
isting MD force- field should describe the total I(k); the 
best agreement was found for the TIP4P-pol2 potential 
from which pair-correlation functions were subsequently 
extracted to represent experiment. However, the inter- 
nal molecular scattering strongly dominates I(k) in x-ray 
scattering and masks the more relevant intermolecular 
scattering such that small, but significant discrepancies 
in phase and amplitude at higher k)^ which determine 
the shape and height of the first 0-0 peak, were not ob- 
served and taken into account. Since the two independent 
studies based on, respectively, neutron and x-ray diffrac- 
tion data arrived simultaneously at similar peak height 
and shape this was understandably taken as proof that 
the 0-0 PCF had been determined correctly; however, 
both studies reproduced in a sense the force-field used 
for the analysis and neither was strictly correct. 

This state of affairs was analyzed more deeply in sub- 
sequent work by Soper, who in two seminal papers^i^ 
first showed that diffraction data do not contain enough 
information to discriminate between structure models of 
strongly different H-bond topology and then that a com- 
bination of x-ray (sensitive to 0-0 and 0-H correlations) 
and neutron diffraction data (sensitive to 0-H and H-H 
correlations) is required to obtain reliable estimates of 
the three PCFs. Considering the significantly reduced 
height of the first 0-0 peak it was concluded that softer 
MD potentials were called for similar conclusions were 
reached based on RMC fits to the same data setsj^iii^ In- 
deed, actually fitting the Hura et al. data set using either 
EPSR^ or RMGiii^ gives a first peak height (2.3) and 
position (2.82-2.85 A) in agreement with the analysis by 
Narten and coworkersi^^iiS^ of their earlier data as well 
as with the Fourier transforms of recent more extended 
data sets,^iHi55 

We now test whether the obtained 0-0 PCF from the 
vdW models, with their low and asymmetric first peak 
at long distance and smeared out second shell, can be 
compatible with the PCF for ambient water under the 
assumption that the interactions are sufficiently well de- 
scribed, but that the simulation protocol may have intro- 
duced too strong constraints on possible structures. That 
HDL-like water would dominate the liquid under ambi- 
ent conditions, i.e. the structure found with the vdW 
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FIG. 6: Mixing of experimental LDL and vdW-DF2 oxygen- 
oxygen PCFs in comparison with PCFs from reverse Monte 
Carlo (RMC) and EPSR analyses (EPSR) of experimental 
data.— 



fuctionals, would be in agreement with what has been 
suggested from x-ray spectroscopy, as well as obtained 
from all scenarios for supercooled waterj ^'^^i-'^^'^^^'^^^ In 
those scenarios fluctuations between HDL and LDL forms 
are assumed and, in view of the vdW functionals seem- 
ingly giving only HDL-like solutions, we explore whether 
adding a "missing" LDL contribution, as postulated in 
these scenarios, could give consistency with experiment 
in terms of the 0-0 PCF. We thus weigh together the 
vdW-DF2 0-0 PCF with a model of LDL to a combined 
PCF and compare with the PCF derived from experi- 
ment using EPSRS and RMCJ^ Since the PBE simu- 
lated structure is far from its preferred density?^ it can 
be assumed to have too large distortions from the "real" 
LDL that could appear as fluctuations in the otherwise 
HDL dominated liquid, and we thus compare with the 
experimental LDL PCF from Soper and Ricci.^^ 

In fitting to the experimental 0-0 PCF we obtain 
agreement (Fig. ^ with a 70:30 mixture between vdW- 
DF2 HDL PCF and the experimentally derived LDL 
PCF.'*^ This ratio is most interesting, since it is very 
close to the original estimate of Wernet et al^ and the 
estimation based on x-ray emission spectroscopyfi2ji£ as 
well as to that from interpreting infrared data in con- 
nection with analysis of a fractional Stokes-Einstein rela- 
tion in water Note furthermore, that quantum effects 
have not been included in the simulations which would 
be expected to bring down and broaden the first 0-0 
correlation additionally ] 

As has been pointed out by Soper-4 when combining 
two separate PCFs one must also consider whether the 
combination introduces additional cross-terms between 
the two, i.e. that the contribution to the total PCF from 
considering pairs of atoms, one from each distribution, 
could change the picture. This would be expected from 
a combination of two highly structured PCFs with well- 
defined peaks occurring at different interparticle separa- 
tions in the two distributions. However, considering that 
both the LDL and HDL local structures give a peak in 
the region of 2.7 - 3 A and beyond that the HDL-like 



PCF is basically without structure it seems likely that in 
this particular case no extra features should be expected 
from cross contributions to a combined PCF. 

The question is naturally why the vdW simulation only 
shows the appearance of HDL-like water and why, in or- 
der to obtain agreement with x-ray diffraction experi- 
ments, it is necessary to artificially add an LDL com- 
ponent. The fact that a combination of an experimen- 
tal LDL 0-0 PCF and that from vdW quite accurately 
reproduces the latest 0-0 PCF of ambient water is of 
course no proof that real water is a combination of the 
two. However, the increased accuracy of the interac- 
tion potential obtained with these latest generation vdW 
functionals indicates that other causes than the non-local 
interaction should be explored to account for the discrep- 
ancy between simulated and measured PCF. 

One potential explanation could be related to the fact 
that the simulation is performed in the NVE ensemble, 
which keeps the volume fixed and thus does not allow 
fluctuations of the density of the box and that this pe- 
nalizes LDL to a greater extent than HDL, once the more 
isotropic vdW interactions are included; the NVE ensem- 
ble is equivalent to adding a pressure to maintain the box 
size which would disfavor fluctuations to lower density 
assuming that the density at ambient conditions corre- 
sponds more closely to that of HDL. The box is further- 
more rather limited with only 64 molecules. In order for 
spatially separated fluctuations between HDL and LDL 
to develop fully it might be necessary to use much larger 
simulation boxes, in particular if the fluctuations are of a 
mean length scale around 1 nm as suggested in Refs. [l^ 
and [3^. There is furthermore some experimental evidence 
from thin water films on slightly hydrophobic surfaces 
that only an HDL related structure is observed even in 
the supercooled regime^^ indicating that if the system 
size becomes very small, indeed only one class of local 
structure is observed and the formation of LDL-like local 
regions is suppressed. 



V. CONCLUSIONS 

The new van der Waals density functionals optPBE- 
vdW and vdW-DF2 show great promise in describing the 
basic structural constituents of liquid water, as seen from 
comparing calculations of water dimer and hexamers 
with benchmark coupled cluster CCSD(T) results i^iS^i^ 
A softening of the structure of liquid water at ambient 
conditions is observed when including vdW interactions, 
consistent with previous work.^^'^^'^^ This is seen from 
the broader angular distributions, the more disordered 
tetrahedrality and asphericity distributions, and from 
the much lower and broader first peak of the oxygen- 
oxygen PCF obtained from the optPBE-vdW and vdW- 
DF2 models compared to PBE. The lower first peak of 
the 0-0 PCF improves the agreement with experiment 
significantly. However, the outer structure is washed out 
by the vdW forces. This has been suggested^ to be re- 
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lated to non-local correlations, but our study of function- 
als with different non-local correlation strength did not 
show any significant difference in the liquid structures, 
while both were found to be very accurate for the water 
dimer. Instead we find that the inclusion of the more 
isotropic vdW interaction shifts the balance over from 
directional H-bonding towards a more close-packed sys- 
tem, i.e. a competition between directional and isotropic 
interactions. 

The vdW simulations seem to be potentially consis- 
tent with a picture of fluctuations between two different 
water structures instantaneously coexisting in nanoscale 
patches albeit not directly observing fluctuations except 
in the sense of obtaining two alternative endpoints with 
vdW forces included (HDL) or excluded (LDL). The rel- 
atively small simulation can only give a picture of the 
local structure of water, and while PBE predominantly 
describes an approximation to low-density water, both 
optPBE-vdW and vdW-DF2, as well as vdW-DF*\ de- 
scribe an approximation to high-density water. By com- 
paring the 0-0 PCFs of the vdW models with PCFs 
from x-rayii^ and neutron^^ diffraction of water at dif- 
ferent pressures we note a resemblance between the vdW 
models and high-density water in terms of effects on the 
second- and third-neighbor correlations while the expan- 
sion of the flrst coordination sphere found in the simula- 
tions may in experiments be counteracted by the pressure 
applied to experimentally generate pure HDL. The com- 
parison to HDL is further supported by the reduction 
of the average number of H-bonds per molecule in the 
vdW MD simulations, which is a result of the isotropic 
vdW forces competing with the directional H-bond for- 
mation. Varying the strength of the exchange interaction 
does not result in a signiflcant change in number of bonds 
once the vdW interaction is included. A 70:30 mixture of 
vdW-DF2 and the experimentally determined LDL PCF 
is compatible with the latest x-ray 0-0 PCF which, how- 
ever, does not constitute proof of a fluctuating real wa- 
ter structure, but indicates the possibility that averag- 
ing over a trajectory obeying less restrictive simulation 
conditions in terms of box size, length of trajectory etc. 
could result in an 0-0 PCF directly comparable with 
experiment. 



Quantum effects are not included in the current simula- 
tions but including them should not qualitatively change 
the consistency with the presented picture. The inter- 
nal Ro_ff bond distance is kept fixed during the simula- 
tions which might affect the hydrogen bonding, but not 
the comparison between PBE and the vdW functionals. 
Lastly, the possibility that the vdW interaction is not 
completely accounted for by the current vdW function- 
als still exists although calibration against various bench- 
marks indicate a quite reliable representation. 

The present work does not resolve the debate on wa- 
ter structure but it suggests for further investigation the 
van der Waals interaction as a physically sound mech- 
anism which affects the balance between directional H- 
bonding and higher packing and may thus indicate a way 
to reconcile the interpretation of recent x-ray spectro- 
scopic data with structures obtained from AIMD sim- 
ulations of liquid water. It is likely that much larger 
and longer simulations in the NPT ensemble are needed 
to determine whether current vdW models support a 
temperature-dependent balance of fluctuations between 
HDL and LDL-like structures in ambient water, as sug- 
gested by recent x-ray spectroscopic and diffraction re- 
sults, and which would be enhanced upon cooling, as 
they must according to all scenarios for water at super- 
cooled temperatures. From the present work it is, how- 
ever, clear that a consistent description of the vdW in- 
teraction in AIMD simulations may possibly provide the 
key to tuning such a balance. 
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